﻿using System;
using System.Numerics;

namespace BassTutor.Analysis
{
    public class FFT : IDFT
    {
        private Reverse rev = new Reverse();

        private void Swap<T>(ref T a, ref T b)
        {
            T x = b;
            b = a;
            a = x;
        }


        public void PermuteInRevBitsOrder(int len, Complex[] data, int bits)
        {
            len /= 2;
            for (int i = 0; i < len; ++i)
                Swap(ref data[i], ref data[rev[i,bits]]);

        }

        private int GetBits(int N)
        {
            int r = 0;
            while (N > 0)
            {
                N >>= 1;
                ++r;
            }
            return r-1;
        }


        public double[] dft(double[] data)
        {
            int N = data.Length;
            int bits = GetBits(N);
            Complex omega, omega_1;
            Complex[] A = new Complex[N];

            for (int i = 0; i < N; ++i)
                A[i] = data[i];

            PermuteInRevBitsOrder(N, A, bits);

            for (int level = 0, group_size = 2; level < bits; ++level, group_size *= 2)
            {
                omega_1 = Complex.FromPolarCoordinates(1, 2*Math.PI / group_size);
                for (int group_begin = 0; group_begin < N; group_begin += group_size)
                {
                    omega = 1;
                    int group_half = group_size / 2;
                    for (int in_group = 0; in_group < group_half; ++in_group)
                    {
                        int near = group_begin + in_group;
                        int far = near + group_half;
                        Complex t = omega * A[far];
                        Complex u = A[near];

                        A[near] = u + t;
                        A[far] = u - t;
                        omega *= omega_1;
                    }
                }
            }

            double[] output = new double[N];
            for (int i = 0; i < N; ++i)
                output[i] = A[i].Magnitude;

            return output;
        }
    }
}
